#include <iostream>
#include<stdlib.h>
#include <time.h>
#include <arm_neon.h>
using namespace std;
const int N=32;//矩阵规模
float A[N][N];
clock_t start, finish;
int count=100;//重复次数
void init()
{
for(int i=0;i<N;i++){
    for(int j=0;j<i;j++)
        A[i][j]=0;
    A[i][i]=1.0;
    for(int j=i+1;j<N;j++)
        A[i][j]=rand();
}
for(int k=0;k<N;k++)
    for(int i=k+1;i<N;i++)
       for(int j=0;j<N;j++)
          A[i][j]+=A[k][j];

}
//平凡算法(串行算法)
void ordinary(){
for(int k=0;k<N;k++){
            float t=A[k][k];
            for(int j=k+1;j<N;j++)
                A[k][j]/=t;
            A[k][k]=1.0;
            for(int i=k+1;i<N;i++){
                for(int j=k+1;j<N;j++)
                    A[i][j]=A[i][j]-A[i][j]*A[k][j];
                A[i][k]=0;
            }
        }
}
//对伪代码第4,5行优化算法
void optimization1(){
for (int k = 0; k < N; k++){
        float t= A[k][k];
        //声明一个包含4个单精度浮点数的向量变量，用A[k][k]初始化
        float32x4_t v1 = vmovq_n_f32(t);
        float32x4_t v0;
        int m;
        for (m = k + 1; m <= N - 4; m += 4) {
            // 从(A[k]+i)地址加载连续4个32位浮点数，保存到v0向量
            v0 = vld1q_f32(A[k] + m);
            v0 = vdivq_f32(v0, v1);//v0/v1
            vst1q_f32(A[k] + m, v0);//把结果装载回原地址去
        }
        //剩余不足四个的部分单独计算
        for (m; m < N; m++)
            A[k][m] = A[k][m] / t;
        A[k][k] = 1.0;
        for (int i = k + 1; i < N; i++) {
            for (int j = k + 1; j < N; j++)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            A[i][k] = 0;
        }
    }
}
//对伪代码8,9,10行优化算法
void optimization2(){
for (int k = 0; k < N; k++) {
        float t = A[k][k];
        for (int j = k + 1; j < N; j++)
            A[k][j] = A [k][j] / t;
        A[k][k] = 1.0;
        for (int i = k + 1; i < N; i++) {
            float32x4_t v1 = vmovq_n_f32(A[i][k]);
            float32x4_t v0, v2;
            int j;
            for (j = k + 1; j <= N - 4; j += 4) {
                v2 = vld1q_f32(A[k] + j);
                v0 = vld1q_f32(A[i] + j);
                v2 = vmulq_f32(v1, v2);
                v0 = vsubq_f32(v0, v2);
                vst1q_f32(A[i] + j, v0);
            }
            for (j; j < N; j++)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            A[i][k] = 0;
        }
    }
}
//两个地方均优化算法
void optimization3(){
for(int k=0;k<N;k++){
    float t=A[k][k];
    float32x4_t v1=vmovq_n_f32(t);
    float32x4_t v0;
    int l;
    for(l=k+1;l<=N-4;l+=4){
        v0=vld1q_f32(A[k]+l);
        v0=vdivq_f32(v0,v1);
        vst1q_f32(A[k]+l,v0);
    }
    for(l;l<N;l++)
        A[k][l]/=t;
    A[k][k]=1.0;
    for(int i=k+1;i<N;i++){
        v1=vmovq_n_f32(A[i][k]);
        float32x4_t v2;
        int j;
        for(j=k+1;j<=N-4;j+=4){
            v2=vld1q_f32(A[k]+j);
            v0=vld1q_f32(A[i]+j);
            v2=vmulq_f32(v1,v2);
            v0=vsubq_f32(v0,v2);
            vst1q_f32(A[i]+j,v0);
        }
        for(j;j<N;j++)
            A[i][j]=A[i][j]-A[i][k]*A[k][j];
        A[i][k]=0;
    }
}
}
//在算法3的基础上对齐
void optimization4(){
for(int k=0;k<N;k++){
    float t=A[k][k];
    float32x4_t v1=vmovq_n_f32(t);
    float32x4_t v0;
    int l;
    for(l=k+1;l<N;l++){
        if(l%4==0)break;
        A[k][l]=A[k][l]/t;
    }
    for(l;l<=N-4;l+=4){
        v0=vld1q_f32(A[k]+l);
        v0=vdivq_f32(v0,v1);
        vst1q_f32(A[k]+l,v0);
    }
    for(l;l<N;l++)
        A[k][l]/=t;
    A[k][k]=1.0;
    for(int i=k+1;i<N;i++){
        v1=vmovq_n_f32(A[i][k]);
        float32x4_t v2;
        int j;
        for(j=k+1;j<N;j++){
            if(j%4==0)break;
            A[i][j]=A[i][j]-A[i][k]*A[k][j];
        }
        for(j;j<=N-4;j+=4){
            v2=vld1q_f32(A[k]+j);
            v0=vld1q_f32(A[i]+j);
            v2=vmulq_f32(v1,v2);
            v0=vsubq_f32(v0,v2);
            vst1q_f32(A[i]+j,v0);
        }
        for(j;j<N;j++)
            A[i][j]=A[i][j]-A[i][k]*A[k][j];
        A[i][k]=0;
    }
}
}
int main()
{
    //平凡算法(串行算法)
    init();
    start=clock();
    for(int q=0;q<count;q++){
        ordinary();
    }
    finish=clock();
    cout<<"ordinary time: "<<((double)(finish - start) / CLOCKS_PER_SEC)/count<< endl;

    //对伪代码中第4,5行的未对齐优化算法
    init();
    start=clock();
    for(int q=0;q<count;q++){
        optimization1();
    }
    finish=clock();
    cout<<"optimization1 time: "<<((double)(finish - start) / CLOCKS_PER_SEC)/count<< endl;

    //对伪代码中第8,9,10行的未对齐优化算法
    init();
    start=clock();
    for(int q=0;q<count;q++){
      optimization2();
    }
    finish=clock();
    cout<<"optimization2 time: "<<((double)(finish - start) / CLOCKS_PER_SEC)/count<< endl;

    init();
    start=clock();
    for(int q=0;q<count;q++){
      optimization3();
    }
    finish=clock();
    cout<<"optimization3 time: "<<((double)(finish - start) / CLOCKS_PER_SEC)/count<< endl;

    init();
    start=clock();
    for(int q=0;q<count;q++){
      optimization4();
    }
    finish=clock();
    cout<<"optimization4 time: "<<((double)(finish - start) / CLOCKS_PER_SEC)/count<< endl;

    return 0;
}
